Phase behavior of a system of particles with core collapse 
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The pressure-temperature phase diagram of a one-component system, with particles interacting 
through a spherically symmetric pair potential in two dimensions is studied. The interaction con- 
sists of a hard core plus an additional repulsion at low energies. It is shown that at zero tempera- 
ture, instead of the expected isostructural transition due to core collapse occurring when increasing 
pressure, the system passes through a series of ground states that are not triangular lattices. In 
particular, and depending on parameters, structures with squares, chains, hexagons and even qua- 
sicrystalline ground states are found. At finite temperatures the solid-fluid coexistence line presents 
a zone with negative slope (which implies melting with decreasing in volume) and the fluid phase 
has a temperature of maximum density, similar to that in water. 



I. INTRODUCTION 

Determination of the phase structure of real materials 
from first principles calculations has been one of the aims 
of statistical mechanics since long ago. Although a qual- 
itative understanding of the processes leading to the dif- 
ferent kinds of phase transitions (between gas, liquid, and 
one or more solid phases) in the pressure-temperature (P- 
T) phase diagram of a classical system has been gained, it 
is clear that the quantitative fitting of the behavior of real 
materials requires a detailed knowledge of the interaction 
between particles and a great deal of computational work, 
that only in recent years has become feasible. 

In addition to the usual materials in which atoms or 
molecules are the basic constituents, in the last years col- 
loidal dispersions have provided a new kind of systems in 
which parameters such as particle size and interaction 
potential can be varied greatlyJd These systems consist 
of a set of latex spheres in colloidal suspension, with 
the aggregate of some amount of non-adsorbing polymer, 
which modifies the interaction potential between the par- 
ticles. Their study has practical importance in relation to 
the properties of many common substances (such as ink, 
paints, cosmetics, blood, etc.). It is clear that a knowl- 
edge of the phase behavior of different model systems is 
important in order to compare the theoretical predictions 
with the experimental results. 

Much effort has been spent in the elucidation of the 
properties of binary mixtures of particles of two differ- 
ent sizes, where segregation, flocculation, gartial crystal- 
lization and other phenomena may occur .□ On the other 
hand, other studies have been directed towards the de- 
termination of the phase behavior of identical particles 
interacting through different model potentials. In this 
case the possibilities for the behavior of the system are 
not so wide as in the case of binary mixtures but, how- 
ever, interesting phenomena occur. It was shown for in- 
stance, that the usual solid-liquid-gas phase diagram of 



particles interacting through a hard core repulsion plus a 
long range attraction isjnodified when the range of the 
attraction is decreased.El More precisely, the liquid-gas 
coexistence curve disappears if the range of the attrac- 
tive potential is lower than about 30% of the hard core 
radius. More interestingly, when the range of the attrac- 
tive potential is reduced below about 8% of the repulsive 
range, a coexistence curve separating two isostructural 
solid phases appear. 

A more obvious isostructural transition occurs in the 
case in which the attractive well is replaced by a repul- 
sive shoulder. In this case, for low pressures, the re- 
pulsive shoulder can sustain a compact structure with 
a lattice parameter related to its range. But when ap- 
plying enough pressure the system must collapse to a 
new compact structure with a lattice parameter given by 
the real hard core of the particles. This kind of mod- 
els, whether with a square shoulder, or a linear ramp 
soft core (which is the one discussed in this paper) have 
been studied since long time ago with the picture of core 
collapse in mindJll Extensions to more general potential 
were also performed.u In recent papers the problem has 
been revisited, and in particular the isostructural transi- 
tion has been studied numericallyJS and analytical results 
have shown that in three dimensions, the ground state 
of a system with a hard core plus a repulsive shoulder 
can be one of various crystalline structures depending on 
parameters!] 

In this paper I show for the hard core plus linear ramp 
model in two dimensions that even the stable zero tem- 
perature structures may be very different from the ex- 
pected triangular structures. The most stable configu- 
ration may be one of a variety of crystalline structures, 
and even a quasicrystal. These structures melt when in- 
creasing temperature. The solid-fluid border in the P-T 
diagram has a zone with negative slope, which implies 
a melting with decreasing in volume and in this region 
the fluid has an anomalous thermal expansion,^ up to a 
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temperature at which a density maximum is attained. 

The paper is organized as foUows. In the next section 
the model is introduced and details of the simulation pro- 
cedure to be used in Section IV are provided. In Section 
III the ground state configurations are analyzed. In sec- 
tion IV, I present detailed results for the P-T phase di- 
agram for a particular value of the parameter a, which 
is defined below. In Section V, possible relevance to real 
systems are discussed and a summary of the results is 
given. 

II. MODEL AND NUMERICAL TECHNIQUE 

The model interaction U (r) between particles that will 
be used here consists of a hard core repulsion at a radius 
To {U (r) |r<ro = zcro for distances larger than a 

value ri, and has a soft repulsive part for < r < ri of 
the form U (r) — Eq (ri — r) / (ri — rg) (Figure |^). This 
interaction gives a model that is a candidate to have an 
isostructural transition between compact configurations 
of lattice parameter tq and ri. Two particles interact- 
ing through this potential in the presence of an external 
force / trying to bring them together, will have a jump in 
the interparticle distance from ri to tq when / exceeds 
the critical value Eq/ (''i — ''o) • This model is preferred 
for numerical simulations instead of the square shoul- 
der model, because it has much less metastabilities when 
varying pressure or temperature. If temperature is mea- 
sured in units of the energy at contact Eq (Boltzmann 
constant is taken to be 1), and distances in units of the 
hard core distance tq, then a = ri/r^ is the only free 
parameter of the interaction potential. 
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FIG. 1. The pair potential used throughout the paper. 
There is a hard core at distance ro and a soft core (linear 
ramp) between ro and r\. Interaction vanishes beyond r\. 

Detailed numerical simulations were performed for a 
two-dimensional system of 256 particles in the NPT en- 
semble, using the Montecarlo-Metropolis technique. A 
trial movement of a particle consists of a displacement 
to a new position chosen randomly inside a cube of a 
linear size of 1% of the mean distance between particles. 



The new position is accepted with a Metropolis algo- 
rithm, considering the energy change due to the move- 
ment. Once each five Montecarlo swepts through all par- 
ticles, a trial global rescaling of all particle coordinates 
and system size is proposed. The rescaling is given by a 
factor chosen randomly within the interval ±0.2%, and 
is done independently for x and y coordinates, in order 
to allow the system to accommodate the different crys- 
talline structures that may appear. A maximum aspect 
ratio for the system of 1.05 is imposed. If the trial vol- 
ume change does not produce hard core particle over- 
lapping, then it is accepted according to the Metropolis 
rule with the value of the energy change A_E, given by 
A£; = P/W - (NT/V) AV + dE. Here N is the number 
of particles, V the volume of the system, dE is the energy 
change associated to the change of interparticle distances, 
and the term — {NT/V) ISV —which assures the correct 
limiting equation of state in the case of an ideal gas (i.e., 
when dE — 0)— accounts for the kinetic energy term that 
—as usual— can be integrated out in the expression for 
partition function and thermodynamic potentials. Dif- 
ferent runs were performed at constant pressure starting 
from random configurations at high temperature, cooling 
down to zero temperature and then warming up. Around 
10000 Montecarlo steps were used for thermalization at 
each temperature, and then 50000 steps were used to cal- 
culate thermodynamic quantities, such as the mean vol- 
ume V per particle at each temperature, the enthalpy per 
particle h, the diffraction pattern of the geometrical con- 
figurations, and the diffusion coefficient of the particles 
D. 

The diffusion coefficient is calculated in the following 
way. The distance traveled by each particle, starting 
from its initial position, as a function of time is recorded, 
and D is taken to be the slope of this function at long 
times. From this definition and the kind of simulations 
performed, it is clear that D tends to a constant at 
high temperatures. The physical diffusion cefficicnt is 
obtained multiplying by temperature. In addition, from 
the diffraction patterns an orientational order parameter 
Bm will also be used. It is defined as 

B„, = J K{k)P{k,9)cx-p{im9)d^k. (1) 

Here P {k, 0) is the intensity of the diffraction pattern in 
the k plane, in polar coordinates, m is an integer chosen 
according the orientational order we are looking for, and 
K [k) is a kernel that cuts off the integral at large k. The 
results are qualitatively insensitive to the form of K (fc) . 
The expression used was K (k) — exp(— fc^). 

Some comments are in order at this point. The use of 
the hard core plus linear ramp potential is motivated —as 
stated before— by numerical reasons. Both the analytical 
results of the following section and the numerical ones of 
section HI do not change qualitatively if a square shoul- 
der, or a parabolic one (with negative second derivative) 
is used instead of the linear ramp. More precise condi- 
tions on the potential are discussed in Section V. Results 
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are presented for the two-dimensional case to clarify the 
discussions of the structures that will be presented in the 
next section. However, the basic properties of the sys- 
tem remain the same in three dimensions. The nature 
of the melting transition in two dimensions is a contro- 
versial point in the literature. However, for rather small 
systems as the one studied here, indications of discon- 
tinuous melting transitions are clearly observed, but is 
not obvious whether they remain when the system size 
goes to infinity. I will speak throughout the paper of con- 
tinuous and discontinuous melting transitions when the 
simulations indicate each case for the particular size used 
in the simulations. 



III. ZERO TEMPERATURE BEHAVIOR 




FIG. 2. Ground state configurations of the system as a 
function of P (= Pro/eo) and a. At the point marked by a 
star, the ground state of the system is a random quasicrystal. 
The black dots in the configurations represent the hard core 
of the particles. 

Since the first numerical works of Alder and 
WainwrigtliEl in the sixties, it has been known that the 
existence of a high enough external pressure P is suffi- 
cient to make a system of otherwise repulsive particles 
to freeze. The minimum energy configuration of a sys- 
tem of particles interacting in two dimensions through a 
potential of the form ~ is a triangular lattice for all 
positive values of 7. It is not so widely recognized the 



great variety of ground states that can be obtained for 
more general potentials, even keeping the restriction of a 
monotonically decaying potential. We will concentrate in 
the already introduced hard core plus linear ramp poten- 
tial (Fig. nl). The ground state configuration of a system 
of particles in two dimensions interacting through this 
potential depends on the values of P and a, and is not 
necessarily a triangular lattice. 

Depending on pressure, nearest particles tend to be at 
distance tq or ri from each other. Intermediate values 
are not preferred because are not energy minima. The 
origin of complex ground state structures in the system 
is related to the competition between two terms in the 
enthalpy H of the system. One is the usual PV term, 
which tends to minimize the volume, and the other is the 
repulsive energy term, which tends to maximize the in- 
terparticle distance. This produces a sort of frustration, 
because both terms cannot be minimized at the same 
time. The two triangular structures with lattice param- 
eters ro and ri (that will be referred to as structures So 
and ^i) correspond to two ways of reducing the enthalpy 
by minimizing one term whereas maximizing the other. 
These are the best compromises in the case of very low or 
very high pressures. However, when both energy terms 
are comparable, lower energy intermediate solutions can 
be found by arranging the particles with a coordination 
number (number of neighbors at distance rg) intermedi- 
ate between and 6 (which correspond to the structures 
Si and 5*0). 




FIG. 3. Tiling of the structures^ S2, and 52 using the 
two-dimensional Penrose tiles, for P = Pqc and a = Oqc. 
The_proportion of thin to fat tiles is 1:3 in 5*3 and 2:1 in 5*2. 
At Pqc, Qqc, any possible random tiling of the plane produces 
a degenerate ground state. 

In fact, different crystalline configurations can be pro- 
posed, and their enthalpy calculated in order to find the 
most stable one as a function of P (= Pr^/eo) and a. 
The result of this analysis is shown in Fig. 0. This figure 
shows the results up to a value of a for which the interac- 
tion to second neighbors in the most compact structure 
(5o) is still zero. The structures in Fig. || were found 
by inspection, and they are the lowest energy configura- 
tions found within each region, but other (more stable) 
structures may have been missed. Note that some of the 
structures have one particle per unit cell (all particles are 
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in translationally equivalent sites), but in others {S3 anc 
S'4) this number is greater than one. For some values 
of a and as a function of P, there are at least three in- 
termediate structure between the triangular ones 5*0 anc 

There is one point in the P-a diagram that deservef 
further discussion, and this is the one marked by a stai 
in Fig. ||. It corresponds to a value of a = aqc = 14 
2sin(18°) ^ 1.618, and P ee Pqc = l/sin(36°) ^ 1.7013 
At this point structures ^2 and ^3 become energeticallj 
degenerated, but more importantly, many other degener- 
ate structures can be constructed. In fact, structure 5*3 
and structure S2 for this value of a may be considered at 
generated by a tiling of the plana-using the two 2D Pen- 
rose tiles as indicated in figure plllj Particles are locatec 
in the far vertices of the thin tile, and in the nearest ones 
of the fat tile. At the point S^^c the enthalpies per parti- 
cle of the two Penrose tiles coincide. Any tiling, with the 
only restriction imposed by the location of the particles 
in the above mentioned vertices of the tiles (these arf 
usually named soft matching rules and the, kind of tiling 
they generate is known as a random tiling)i2l'E£l) generates 
a possible ground state of the system. The proportion o: 
thin to thick tiles used to construct the ground state if 
arbitrary (as long as the soft matching rules can be sat- 
isfied). For proportions close to the value corresponding 
to a perfect Penrose lattice (nearly 1) the structure ob- 
tained is a random quasicrystal. For pressures lower thar 
Pqc the structure with a maximum fraction of thin tiles 
(^2) is preferred, because the thin tile has lower enthalpy 
On the contrary, for pressures higher than Pqc the pre- 
ferred structure is that with the largest proportion of fal 
tiles (5*3). Quasicrystalline ground states are only stabk 
at the point Sqc- However, the soft matching rules al- 
low for many ways of generating them (compared to the 
more rigid structures S2 and S3), and this implies thai 
at finite temperature the quasicrystalline structure wil 
be stabilized due to entropic effects, as we will see in the 
next section using numerical simulations. 

IV. PRESSURE-TEMPERATURE PHASE 
DIAGRAM 

After having discussed its ground state properties, we 
focus now on the behavior of the system at finite tem- 
peratures. The main interest is in the localization of the 
fluid-solid transition line. I will present now the numer- 
ical results obtained at different values of pressure, ir 
the particular case a = 1.65. The system is initializec 
in a random configuration at high temperature (well in- 
side the fluid phase) and the temperature is progressivelj 
reduced to zero, and then increased again. 
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FIG. 4. Volume v, and enthalpy h (= E/N + Pv) per par- 
ticle, sixfold orientational order parameter Be, and diffusion 
coefficient D of the system as a function of temperature for 
P = 0.5, for a swept decreasing and increasing temperature 
{Be and D are given in arbitrary units, T is in units of eo). 
Note the difference in the snapshots of the system at the 
same temperature on heating and cooling, within the hys- 
teresis loop. 

In the ranges of pressure in which the structures Sq and 
Si are the most stable zero temperature configurations 
(according to Fig. H), a sharp solidification transition is 
obtained when reducing temperature. This can be seen 
in figures and ^ for three different values of pressure 
within this range. In the first two the solidification is 
into the 5i structure, and for the third one into the So 
structure. In the three cases the solidification transition 
is clearly identifiable by the hysteresis loop in the volume 
or the enthalpy of the system, which coincides with the 
vanishing of the diffusion coefficient and the appearance 
of a finite sixfold symmetry of the diffraction pattern. 
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FIG. 5. Same as Figure ^ for P = 



1.0. 



An additional check of the existence of a sharp sohd 
fluid transition may be obtained through a long Simula 
tion at the equilibrium temperature between solid anc 
fluid. In this case, the volume of the system should flue 
tuate between two clearly different values corresponding 
to solid and fluid phases. Results of this simulation ar( 
shown for the case P = 1, and T = 0.082 in Fig. 0. Foi 
this simulation the system was initialized in a fluid equi 
librium configuration at the corresponding values of J 
and T. After about 5x10^ Montecarlo steps the systen 
jumps to the solid phase. After around 2.7x 10^ steps th( 
systems makes a new short transition to the fluid state. 



FIG. 6. Same as Figure g for P = 6.0. 
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FIG. 7. Time evolution of the volume per particle v for 
P = 1.7 and T — 0.082, close to the fluid-solid transition. 
After about 5x 10^ Montecarlo steps the system jumps to the 
solid phase. After around 2.7x10^ steps the systems makes a 
new short transition to the fluid state. 

The most important characteristic to be noted in Figj^, 
is that the melting occurs with a reduction in volume for 
this value of pressure. In addition, the fluid right af- 
ter melting has also anomalous thermal expansion up to 
some temperature at which a density maximum is at- 
tained. These characteristics imply a negative slope of 
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the solid-fluid coexistence curve, which is in fact obtainec 
from the simulations as we will see later. The compressive 
melting of the system in this region has its origin in the 
fact that the usual volume reduction when temperature 
is reduced is overcome by the expansion produced wher 
particles diminish their kinetic energy and move out ol 
the soft core of their neighbors. Illustrating this effect, ir 
Fig. ^ we see snapshots of the system at different temper- 
atures passing through the liquid-to-solid transition. Ir 
the fluid phase there are particles at distances lower thai 
ri from each other, whereas in the solid phase the min- 
imum distance between particles is ro (except for some 
defects in the structure, which appear mainly because 
of the impossibility of accommodate 256 particles in e 
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FIG. 9. Volume per particle v, enthalpy h, and diffusion 
coefficient D of the system as a function of temperature for 
P = 1.3, for a swept decreasing and increasing temperature. 
The zero temperature configuration reached and its diffrac- 
tion pattern are also shown. 



FIG. 8. Snapshots of the system for P — 1.7 when de- 
creasing temperature, through the fluid-solid transition, illus- 
trating the anomalous freezing. In the fluid state there are 
particles at distance lower than n , whereas in the solid phase 
all particles (except a few defects) are located at distance rj 
from their neighbors. The temperature at the second pane' 
corresponds to the maximum density of the system (see Fig 

!)• 



The phases iS'2, 5*3, and S4, expected to be the ground 
states of the system at intermediate pressures for a — 
1.65, are not straightforwardly obtained in the simula- 
tions. Instead, rather disordered states are obtained. In 
figures I, |l^, and^ we can see the magnitudes v, h, and 
D for pressures P = 1.3, 1.7, and 3.8, together with the 
zero temperature configuration found and the diffraction 
pattern of the zero temperature structure. In all the in- 
termediate pressure range (1.2 < P < 4) the enthalpy at 
zero temperature obtained in the numerical simulation 
is never lower than the corresponding to the expected 
ordered structures of Fig. || as it is shown in Fig. |l^, 
indicating that probably the configurations of Fig. ^ are 
really the fundamental states, but they were not reached 
in the simulations. The configurations obtained in the 
simulations reflect the equilibrium states at some finite 
but small temperature, where entropic contributions to 
the free energy are important, and they are metastable 
at zero temperature (note the existence of chains in Fig. 
^, and the pentagons and hexagons in Figs. fo| and 
Looking at the diffraction patterns in Figs. |l^, and 
|ll| , some of them {P — 1.3, and P = 3.8) show no sign 
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of orientational order. Others {P = 1.7) clearly indicate 
a tenfold symmetry, characteristic of a quasicrystal. Ir 
the cases where the low temperature state has no ori- 
entational order, the volume or enthalpy of the systerr 
do not show any abrupt solidification transition. In the 
cases where the low temperature state has rotational or- 
der the volume and enthalpy of the system show a smal 
hysteretic behavior, suggesting an abrupt solidificatior 
transition. To confirm this fact, long runs were per- 
formed at the expected transition temperature, recording 
the temporal evolution of volume and enthalpy. An ex- 
ample of the results for the case P — 1.7 is shown in Fig 
p^ . The histogram shows a clear bimodal distributior 
between two values corresponding precisely to the values 
of enthalpy expected from Fig. 10 at this temperature. 





FIG. 10. Same as Fig. g for P = 1.7. In this case the 
tenfold orientational order parameter Bio is also shown ir 
the last panel. Note the small histeresis loop in v and h. 



FIG. 11. Same as Fig. | for P = 3.8. 

We saw in the previous section that at zero temper- 
ature the quasicrystalline structure is stable only at a 
particular value of P and a. At finite temperatures this 
structure is stabilized due to entropic effects, because 
there are many ways of construct the random tiling, fa- 
voijiag this structure against the more rigid ones 5*2 and 
S'aJlj A related model of a quasicrystal using two kinds 
of particles of differeiilj-S|izes has been studied by Strand- 
burg and coworkers .EiJo In our case, the quasicrystalline 
state is obtained in a system of only one kind of particles. 

For other values of pressure, the smooth solidification 
of the system, and the absence of any obvious order in 
the low temperature structures obtained, suggest that 
the system freezes into a glassy state. However, more 
detailed calculations of the diffusion coefficient and other 
magnitudes in larger systems are needed to confirm this 
point. 
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FIG. 12. Ground state enthalpy per particle as a function 
of P from the simulations (black dots) and the analytical ex- 
pression for the possible ordered structures. All points lie 
above at least one of the lines corresponding to the ordered 
structures. 



in the other cases is also indicates, and in this case the 
error bars indicate the approximate temperature range 
in which the diffusion coefhcient changes between 10 and 
50 % of its value at high temperatures. 
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_ FIG. 13. Time evolution of the enthalpy per particle for 
P = 1.7, T = 0.0505, close to the fluid-quasicrystal tran- 
sition. The histogram shows a clear double peak structure, 
with mean values compatible with those obtained from Fig. 
|lo| for this temperature. 

The numerical results are summarized in the phase di- 
agram of Fig. The sharp fluid-solid transition in the 
case of structures 5*0 and S'l are shown, as well as the cor- 
responding to a quasicrystalline state. The error bars in 
these cases are taken as the width of the hysteresis loop 
in the enthalpy or the volume of the system. In addition, 
the approximate temperature where the system freezes 



FIG. 14. Pressure temperature phase diagram from the 
simulations (dsd stands for disordered). See text for more 
details. 

We know from the results of the previous section that 
this phase diagram (particularly at intermediate pres- 
sures) is metastable at low temperatures. Since the nu- 
merical simulations are not able to reach the fundamental 
state in some cases, the determination of the equilibrium 
phase diagram at all values of P and T requires the direct 
comparison of the free energies of the structures found 
in the simulations and those known to be more stable in 
some cases (structures ^2, 5*3, and 5*4, for a = 1.65). The 
Gibbs free energy of the system can be obtained from the 
numerical simulations through the formula 

I'LdP-^dT (2) 
T2 Ti T r2 ^ ^ 

where 1 and 2 stand for two set of values Pi, Ti and 
P2, T2, and the integration is through an arbitrary (re- 
versible) path in the P-T plane joining points 1 and 2. 
This determines the free energy of those structures ob- 
tained in the numerical simulations up to an overall con- 
stant. Also the free energy of structures 5*2, S'3, and 
Si may be determined in this way, by setting up the 
configuration of the system at zero temperature in these 
structures, and then performing a numerical simulation 
increasing temperature. After that, all that remains to 
be able to compare the free energies is to fix the additive 
constants. This was done by introducing in the model an 
additional external potential characterized by a strength 
W , with a periodicity chosen to favor the formation of 
the required structure. The reversible path from the or- 
dered structure to that obtained in the simulation for a 
given point Pq, Po consisted of four steps: increasing W 
from zero to some large value, increasing T from Pq to a 
large value, decreasing W down to zero, and decreasing 
P down to Pq. The difference in Gibbs free energy AG 
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between the two structures was calculated through this 
path by a generalization of the formula , given by 



A(G/T)|p„,To 
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(3) 



where Ew is the potential energy of the particles in the 
artificial external potential, and f indicates the integra- 
tion along the above mentioned path. This was done 
three times with different external potential to fix all ar- 
bitrary constants between free energies of structures 5*2, 
53, S'4, and the ones obtained in the simulations. After 
that, the free energies can be compared and the thermo- 
dynamical phase diagram constructed. 

The complete thermodynamic phase diagram is shown 
in Fig. [ij, where the stability region of each phase is 
shown. It is seen that the quasicrystalline state is in fact 
thermodynamically stable in a finite range of P and T, in 
spite of the value of a (= 1.65), which is not the optimum 
value (ckqc = 1.618) for the quasicrystalline structure. 
Only in the case a = ofqc the quasicrystal is stable down 
to zero temperature, at the point P = P^c = 1.7013. 
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FIG. 15. Complete pressure-temperature phase diagram 
for a = 1.65. The errors in the limit of structures 5*2, S3, 
and 54 are estimated to be ±0.01 in temperature. See text 
for discussions. 



V. SUMMARY AND DISCUSSION 

In this paper I have discussed the phase behavior of a 
classical model of particles interacting through a particu- 
larly chosen isotropic potential. The interaction consists 
of a hard core plus a linear ramp potential, producing 
an effective greater size of the particles at low energies. 
The ground state of the system shows different periodic 
arrangements of the particles, depending on the values 
of P and a. There is a particular value of P and a at 
which the ground state is a random quasicrystal. These 
periodic structures melt when temperature is increased 
through discontinuous phase transitions. In addition, at 



finite temperature, the quasicrystalline state is stabilized 
due to entropic effects. 

It is usually believed that "...two-dimensional 
monoatomic systems interacting-, with central forces al- 
ways form a triangular lattice. "t3 Although this is so for 
power laws and other kinds of interactions, the ground 
state configurations of our system show that this is not 
true in general. Even considering only crystalline ground 
states, there may be more than one atom per unit cell 
(2 for S3 and 5 for S'4) and unequivalent sites within the 
structure (2 for ^4). 

Having established the properties of the system for our 
model interaction, it is of basic importance to identify 
the conditions that an interaction potential must satisfy 
in order to obtain the kind of structures we got in our 
model. Although it is rather difficult to solve the problem 
in general, something can be said about. It is clear that 
the non-analyticity of the potential used is not a crucial 
element, actually, analytic potentials arbitrary close to 
the one used here may be easily constructed. First of 
all, let us analyze the necessary condition for the local 
stability of a triangular structure. Consider particles in- 
teracting through a potential U (r) , and distributed in 
a triangular lattice of lattice parameter a. At any lat- 
tice site the potential created by all other particles must 
have a minimum in order for the structure to be (locally) 
stable. The potential on that site (taken to be r = 0) cre- 
ated by a particle at a generic position r is, up to second 
order, of the form U" (r) + ^^y^ + U' (r) x + U (r), 
where x (y) is the coordinate along (perpendicular to) 
r. This potential must be summed up for all particles, 
and for lattices with rotational symmetry C3 or higher it 
must reduce to a isotropic form. Considering the invari- 
ance of the trace of cuadratic forms under rotations, the 
cuadratic part U2 (r) of the final effective potential can 
be written as 



t^2(r)=^^ 



U" id. 



U' id,) 



(4) 



where the sum is over all other particles, located at dis- 
tances di (c?i — a), and Ui is the number of particles at 
those distances (for three dimensions the term containing 
U' gets an additional factor 2). The positiveness of this 
form is the condition for the stability of the lattice under 
small displacements of a single particle (global stability is 
more difficult to characterize), supposing the lattice pa- 
rameter is fixed from outside. For U (r) l/r'^ all terms 
of the sum in (^) are positive, and the structure is locally 
stable. For our potential, and analyzing a structure with 
lattice parameter a, ri > a > tq , nearest neighbors have 
U" (r) = and U' (r) < 0, and the structure is unstable. 
If the stability condition (^ is not satisfied for some lat- 
tice parameter a, then it indicates at least the existence 
of an isostructural transition as a function of pressure be- 
tween two triangular structures with lattice parameters 
ao <a and oi > a. However, if first neighbors interaction 
dominate, it is easy to see that a S2 kind of structure is 
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more stable. In fact, consider two triangular structures 
with lattice parameters qq and ai which are degenerate 
at some pressure P. Their enthalpy per particle must be 
equal. This implies —if only first neighbors interaction is 
important— that 

PalVi/2 + 3U{ao) = PalV3/2 + 3U{ai). (5) 

Now, the enthalpy of a 5*2 structure with nearest neigh- 
bors distance oq, and second neighbors distance ai is 
given by 

Pao^aj - all A + {/(ao) + 2L/(ai), (6) 

and it is easy to see that this number is lower than the 
corresponding to the triangular structures. We conclude 
that stable structures other than the triangular one will 
occur if the condition 

U" (r) + < (7) 

r 

is satisfied for some value of r. In cases where next- 
neighbors interactions cannot be neglected, a case by case 
analysis based on Eq. (Q) is needed. 

Our model does not contain an attractive part in the 
potential, and for this reason it lacks a liquid phase. How- 
ever a liquid phase can he obtained within a generalized 
van der Waals approachJl3 in which an attractive interac- 
tion is included through an energy term proportional to 
—v~^. In addition to the appearance of a liquid-gas coex- 
istence line, this modification only renormalizes pressure, 
and does not affect the structure of the phases, nor the 
nature of the anomalies of the phase diagram that were 
discussed above. 

The qualitative features of the phase diagram in the 
two-dimensional case have also been obtained in simu- 
lations with a three-dimensional system. The minimum 
value of a necessary to get the volume anomaly at melt- 
ing is about 1.2, both for two- and three-dimensional sys- 
tems. 

The results presented here might have importance in 
another context. The anomaly in the fluid-solid coexis- 
tence line, and the negative thermal expansion coefficient 
in this region remind strongly the same effects occurring 
in water. Actually, an interaction potential with a dou- 
ble minimum (which is related to the potential consid- 
ered here) in one dimension has been analyzed recently, 
and has been suggjested to be the origin of the density 
anomaly in water ,1111 but this proposal has been ques- 
tioned since it pjsms to give anomalous behaviors only in 
one dimension.113 The results presented here show clearly 
that simple models with anomalous behavior may be con- 
structed in two and three dimensions. A comparison of 
the phase diagrams of water with the corresponding to 
our system reveals in fact many coincidences, as the men- 
tioned anomalies and the position of the different crys- 
talline phases in the P-T diagram, which occur near the 
pressure where the melting temperature is minimum. In 



addition, amorphous structures in ice are well known,t3 
and the underlying mechanisms responsible for their for- 
mation may be related to those that originate our dis- 
ordered structures when cooling down at intermediate 
pressures. 

Finally, in order to get a better understanding of the 
dynamic and thermodynamic properties of the kind of 
system we dealt with, it would be interesting to find an 
experimental realization of an interaction potential sat- 
isfying condition (^. A colloidal system seems to be the 
natural candidate where to look for that realization. 
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